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We studied the shock propagation and its stability with the causal dissipative hydrodynamics in 
1+1 dimensional systems. We show that the presence of the usual viscosity is not enough to stabilize 
the solution. This problem is solved by introducing an additional viscosity which is related to the 
coarse-graining scale of the theory. 
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I. INTRODUCTION 

It is by now widely accepted that the basic features of 
collective motions in relativistic heavy-ion collisions can 
be well described as those of a hydrodynamical motion of 
an (almost) ideal fluid Several studies on the effects 
of viscosity are available in literatures and support such 
a vision @, H, 0, @ ■ In addition, the possibility of the 
existence of the lower bound of the shear viscosity co- 
efficient has been discussed by assuming the AdS/CFT 
correspondence @. 

However, strictly speaking there are still several open 
questions in hydrodynamic approaches of the heavy-ion 
collisions. Hydrodynamic observables are not so restric- 
tive to determine uniquely many unknown factors, such 
as initial condition, equation of state, dissipation mech- 
anisms, and so on. In particular, even the theory of a 
relativistic dissipative fluid is not yet well understood. 
Therefore, in order to conclude that the matter created 
in relativistic heavy-ion collisions really behaves just as 
an ideal fluid, we need to investigate the effect of dissi- 
pation more carefully. 

The difficulty of the construction of the relativistic 
dissipative hydrodynamics is because of the problem of 
acausality and instability: a naive relativistic general- 
ization of the Navier-Stokes equation present an infinite 
propagation speed of pulse signals and the solution is un- 
stable for small perturbations. To solve this problem, one 
should, for example, take into account memory effects by 
introducing a relaxation time [![. 

So far, there are several attempts to study the effect 
of dissipation to relativistic heavy-ion collisions by im- 
plementing the 1+1 and 2+1 dimensional calculations 
4]. In these studies, they mainly deal with cases where 
the deviation of the ideal fluid is rather small. On the 
other hand, there is an evidence that the bulk viscos- 
ity becomes large in the vicinity of the critical point Q. 
In addition, to know the limitation of the theory itself, 
we should investigate the behavior of solutions for large 
viscosity as well. In particular, for higher energies such 
as in the LHC regime, we expect that the relative im- 
portance of the viscosity becomes significant. This point 
is essential since the problem of acausality or instabili- 
ties is directly related with the size of viscosity and the 
relaxation time mentioned above. 

Other interesting aspects of viscous fluid dynamics ap- 



pear when discontinuities emerge during the time evolu- 
tion or already exist in the initial condition. In the usual 
application of hydrodynamics, only very smooth initial 
distributions, both in energy and velocity, have been ap- 
plied. In such cases, the fluid profile usually remains 
smooth in time and no special attention is required for 
the treatment of discontinuities. 

However, it is sometimes necessary to discuss the ex- 
treme cases which involve discontinuities. For example, 
the Landau type initial condition is often discussed as an 
interesting possibility of meson production mechanism 
with null initial velocity in the p-p collisions and A-A 
collisions [8(. Furthermore, the shock phenomena, which 
are typical discontinuous propagation of hydrodynamical 
variables, also may occur in the heavy-ion collisions by 
high energy jet propagations in the QGP 0, @- We also 
expect a formation of shock wave in the region near the 
coexisting phase, if the QCD-hadron phase transition is 
of first order, where the velocity of sound vanishes (or 
becomes very small). Shock phenomena in relativistic 
heavy ion collisions, if any, are particularly interesting 
since they would furnish genuine hydrodynamical signals. 

It is well-known that to deal with dynamical discon- 
tinuities, such as shocks, is not a simple problem. We 
have to introduce some specific techniques such as the 
Godunov method or artificial viscosity (pseudo- viscosity) 
[To| to achieve physically meaningful results interpolating 
smoothly the discontinuity. For a causal dissipative hy- 
drodynamics, there are a few works where the dynamics 
of shock discontinuities are discussed [Til. However, 
to the author's knowledge, detailed numerical study of a 
causal dissipative fluid dynamics involving discontinuities 
has not yet been done. 

In this paper, we study in detail the dynamics of a vis- 
cous fluid. For the sake of simplicity, we concentrate our- 
selves to 1+1 dimensional systems. The basic objective 
is to analyze the problem of instabilities associated with 
causality through several numerical examples, in partic- 
ular, which contain discontinuities or shock phenomena. 
We also pay attention to the origin of the so-called "ar- 
tificial viscosity" for the numerical calculations of shock 
phenomena in the framework of causal viscous hydrody- 
namics. We investigate its role and interpret its origin in 
terms of the scale of the coarse graining introduced in the 
hydrodynamical theory. To clarify this point, we use the 
smoothed particle hydrodynamics (SPH) as the numer- 
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ical method instead of the commonly used space-fixed 
grid methods. Furthermore, having in mind the collec- 
tive QGP motion in the LHC energy regime, we restrict 
ourselves to systems described by an equation of state of 
a baryon-number free, massless relativistic gas. 

In the next section, we briefly review the causal dis- 
sipative hydrodynamics and describe how the memory 
function is introduced. To check causality and stability 
of our theory, we discuss the dispersion relation for the 
propagation of a perturbative plane-wave. In Sec. Ill, we 
describe the SPH method applied to our problem which 
introduces the coarse-graining scale h. In Sec. IV, we 
show several examples to reveal the effects of viscosity, 
in particular, the case of steady shock wave propagation 
induced by a large pressure gradient in the initial condi- 
tion. In the first example, we show the universal relation 
between pressure and viscosity in the fluid expansion to 
vacuum. This relation is satisfied independently of ini- 
tial conditions and equation of states. In these examples, 
some quick oscillating modes appear in the dynamics in- 
volving shocks, leading to instabilities of the numerical 
solutions. That is, the normal viscosity does not neces- 
sarily damp all the high frequency modes. The presence 
of such modes indicates the necessity of a new ingredient 
for the theory to be physically consistent. In Sec. V, we 
introduce an additional viscosity, having a different scale 
of the memory function associated with the coarse grain- 
ing size, which solves this problem. We show several 
examples where the dynamics of shock phenomena are 
described satisfactorily in this scheme. Finally in Sec. VI, 
we summarize our work on the analysis of a causal dissi- 
pative hydrodynamics in 1+1 and discuss the problems 
still open in such theories. 



II. RELATIVISTIC DISSIPATIVE 
HYDRODYNAMICS 

Various theories have been proposed to incorporate dis- 
sipation consistent with causality and stability; the di- 
vergence type theory [TI|, the Israel-Stewart theory [l4| 
and its extension based on the extended irreversible ther- 
modynamics Jk| i Carter's theory [l6j], Ottinger-Grmela 
formulation [l7|| and the memory function method 

Here we briefly review the memory function method 
[H to obtain a causal dissipative hydrodynamics for 1+1 
dimensional case. That is, we ignore the motion in the 
transverse direction and concentrate on the longitudinal 
dynamics. 

As was mentioned in the introduction, we just consider 
the case of vanishing baryon chemical potential. In this 
case, the hydrodynamical equation of motion can be writ- 
ten only as the conservation of the energy-momentum 
tensor, 



0. 



(1) 



together with the thermodynamical relations. We use the 
Landau definition for the local rest frame and assume, as 



usual, that the thermodynamic relations are valid in this 
frame. Then the energy-momentum tensor is expressed 
as 



= (e + p + n) - (p + n) gH», 



(2) 



where, e, p, and LT are, respectively, the energy den- 
sity, pressure, four velocity and bulk viscosity. 

In the presence of these irreversible currents, the en- 
tropy is not conserved anymore. Instead, from Eq. (TT]), 
we have UM 



(3) 



where the entropy four-flux is identified by Landau and 
Lifshitz as 



(4) 



In irreversible thermodynamics, it is interpreted that 
entropy production is the sum of the products of thermo- 
dynamic forces and irreversible currents. From Eq. ([3]), 
we define the thermodynamic force as 



F = d a u a . 



(5) 



To satisfy the second law of thermodynamics locally, and 
hence the positiveness of the entropy production, Landau 
proposed that the irreversible current should be propor- 
tional to the thermodynamic force [18l |. 



n = -(F = -(d a u a , 



(0) 



where C is the viscosity coefficient. However, it is known 
that the derive d eq uations have the problem of acausality 
and instability pjj H(| [Hj] . To solve these difficulties, we 
introduce a memory effect to the irreversible current by 
using a memory function One of the simplest forms 
of the memory function is 



G(r,r') 



ft> 



TR (T 



dr" 



(7) 



The relaxation time tr(t) is, in general, a function of the 
local proper time r = r (r, i) through the thermodynam- 
ical quantities. Then, the irreversible current is modified 
as follows; 

n (r) = - [ dr'G (t, t') (d a u a (/) + e-^ T - Ta ^ Tli n Q , 

(8) 

where Ho is the initial value given at tq. 

Because of the modification of the relation between the 
irreversible currents and the thermodynamic forces, the 
algebraic positivity of the second law of thermodynamics 
is not satisfied. However, we checked that the second law 
of thermodynamics is not violated in all examples dis- 
cussed in this paper numerically. It is worth mentioning 
that the algebraic positivity is not automatically satis- 
fied even in the Israel-Stewart theory (See the discussion 
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below Eq. (2.31) in [Tj]). To solve this problem, the con- 
cepts of thermodynamics have to be extended. See [l5[ 
for details. 

The integral expression ([8]) are equivalent to the fol- 
lowing differential equation, 



n 



-(d a u a - t r 



dn 



(9) 



M is the total derivative with respect to 
In the 1+1 dimensional case discussed 



where d/dr = u^d, 
the proper time 
here, the equations derived above are equivalent to those 
of the Israel-Stewart theory. Thus, the conclusions in this 
paper is applicable also to the Israel-Stewart theory. 



A. Propagation speed of signals 

We discuss the propagation speed of the 1+1 dimen- 
sional system. For this purpose, we consider a small per- 
turbation of the three independent variables in the form 
of a plane wave, 



5e 

fijjl I — ^iujt—ikx 

su 



propagating in a hydrostatic equilibrated background. 
Then, the linearized hydrodynamic equation for these 
perturbations should satisfy 

iuj —ik(e + p) 
a(—ik) iu(e+p) 
-ikC 

. . (10) 

where a = dp/ de is square of the velocity of sound in the 
absence of the viscosity. 

To have the non-trivial solution for the perturbation, 
the determinant of the 3x3 matrix in the left-hand side 
of Eq. lfTTj)) should vanish so that ui should satisfy the fol- 
lowing dispersion relation, 




Lo - ak 2 



(e + p) 1 + itruj ' 
whose solution for uj can be written as 

A 1 i 
u = x + + o — ' 

where 



(11) 
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FIG. 1: The real part of the frequency ui as a function of k. 
There are two propagating modes (dashed lines A and B) and 
one non-propagating mode (solid line C). 
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FIG. 2: The imaginary part of the frequency uj as a function 
of k. The two propagating modes A and B shown in Fig. 
[1] are degenerated (dashed lines). The solid line is the non- 
propagating mode C. 



with 



l _ c i 



b r R e + p 



(15) 



The asymptotic forms of the dispersion relation for 
large k are 



±k\/f + n 



2r R (l+ab) 

i M M + O (k- 1 ) 

T R (l+ab) \ ) 



O k- 



(16) 
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whereas for small k, 




y/ a -( a +l/4b)r^/b lk 2 TR 

{l-T%k'/b) + 2{b-rik') . (17) 

i/TR 



We have three solutions; two of them are propagating 
modes and the remaining one is a non-propagating mode 
(u> pure imaginary). From Eqs. (jl6| and (|17| we can see 
that the sound velocity of propagating modes reduces to 
that of the ideal one for small k, on th e other h and, for 
large k, the sound velocity is given by ^/l/b + a. In Figs. 
[T] and [2J we show, respectively, the real and imaginary 
parts of oj from these dispersion relations as functions of 
momentum k for a = 0.1, b = 6, and the temperature T = 
200 MeV. In Fig. [U the lower two curves (solid and dot- 
dashed lines) correspond to the upper line of Eq. l(TC|) and 
the horizontal line (dotted) corresponds to the second line 
of Eq. (|16p . As we see from Fig. [21 the imaginary part of 
d is always positive for all k and converges to constant 
values. Positivity of the imaginary part guarantees the 
stability of the plane-wave perturbation. 

We assume that the propagation of physical quantities 
for the propagating modes are characterized by the group 
velocity. Then, from this relation, the maximum velocity 
in such a theory is determined as 

v M = y/l/b + a. (18) 

It should be noted that Eq. (|T8j) gives the velocity of 
sound in the causal hydrodynamics. As a matter of fact, 
in the vanishing £, this definition coincides with the ve- 
locity of sound of ideal fluid. It is clear that the group 
velocity becomes infinite at the limit of tr — > 0. In such 
a situation, always some portion of the matter tries to 
propagate with velocity larger than the velocity of the 
light for any initial condition. However, since the hy- 
drodynamic equation is covariant, the presence of the 
light-cone singularity forbids such a propagation. Then 
the matter tends to accumulate at the light-cone and in 
such situation the linearized wave analysis in the homo- 
geneous back ground at rest will breakdown. This kind 
of conflict between the causality and relativistic covari- 
ance leads eventually to the instabilities of the solution 
near the light-cone [iS E3, EU . 

As for the non-propagating mode, we find that the 
imaginary part becomes a constant for large k. That is, 
the large k component just damps exponentially. This 
is different from the case of a diffusion process where 
the non-propagating mode behaves as k 2 in the large k 
limit, which leads to an infinite propagation speed. In 
this sense, our theory is causal. 

B. Parameters 

There are so far two approaches to estimate the trans- 
port coefficients; kinetic approach and microscopic ap- 
proach. The calculation of the bulk viscosity coefficients 



are much involved in the kinetic approach. We cannot use 
the Boltzmann equation since it contains only the infor- 
mation of two-body collisions. We have to use the (mod- 
ified) Enskog model or the Bogoliubov-Cho-Uhlenbeck 
equation where the multiple collision effects are included. 
In the microscopic approach, it is known that the trans- 
port coefficients are calculated by using the Green-Kubo- 
Nakano (GKN) formula. However, the GKN formula is 
the formula for the Newtonian fluid like the Navier-Stokes 
fluids, and hence we cannot use for the transport coef- 
ficients of the causal dissipative hydrodynamics because 
it is a non-Newtonian fluid. To calculate the transport 
coefficients of the causal dissipative hydrodynamics, we 
have to derive a new formula. One possibility of such 
a new formula is proposed in [13, ES EH- In any case, 
no reliable estimate for the bulk viscosity coefficient is 
available. 

In the present analysis, we will not deal with a precise 
quantitative description of the behavior of the matter 
created in heavy-ion collisions but rather interested in 
qualitative role of viscosity. Thus, we assume a simple 
expression for the bulk viscosity coefficient £ as usually 
adopted for the shear viscosity, 

C = as, 

where s is the entropy density in the local rest frame. 

Another important parameter of the theory is the re- 
laxation time tr. As was mentioned in the previous sec- 
tion, causality constraints the relation between the two 
parameters, C and tr through vm < 1 Thus, for the sake 
of simplicity, we parametrize the relaxation time by tak- 
ing b as constant. This determines the relaxation time 
as 



The later examples are presented in terms of these pa- 
rameters a and b. However, there is no theoretical reason 
that a and b are constants, but they may depend on ther- 
modynamical quantities (see the later discussion) . In the 
following calculations, we use two values for the param- 
eter a, 0.1 and 1, with fixed value of b = 6. 



III. SMOOTHED PARTICLES 
HYDRODYNAMICS 

To solve numerically the relativistic hydrodynamic 
equations we use the Smoothed Particle Hydrodynamic 
(SPH) method. This method was initially introduced 
[25] for application in astrophysics. Using variational ap- 
proach, this method was extended to the application for 
heavy ion collisions (26|. 

The original idea of the SPH method is to obtain an 
approximate solution of hydrodynamics by parameter- 
izing the fluid into a set of "effective particles". How- 
ever, in the application to the heavy ion dynamics, the 
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SPH method is not a mere mathematical discretization 
scheme, but can be interpreted as a physical model of 
the collective motion in terms of a finite set of dynamical 
variables. 

To see this, let us consider a distribution a (r, t) of any 
extensive physical quantity, A. In a system like the hot 
and dense matter created in heavy ion collisions, the be- 
havior of a(r, t) contains the effects of whole microscopic 
degrees of freedom. We are not interested in the ex- 
tremely short wavelength behavior of a (r, t) but rather in 
global behaviors which are related directly to the exper- 
imental observables. Therefore, we would like to intro- 
duce a coarse-graining procedure for a. To do this, we in- 
troduce the kernel function W (r — f , h) which maps the 
original distribution a to a coarse-grained version ace 
as, 

a CG (r, t) = J a(r, t)W (r - f , h) dr (19) 



where W is normalized, 



J W(r,h)dt = 1, 



(20) 



and has a bounded support of the scale of h, 



satisfying 



W(r,h)^0, \r\>h, 



limW(r.fc) =5(r). 



Here, h is a typical length scale for the coarse-graining 
in the sense that the kernel function W introduces a cut- 
off in short wavelength of the order of h. Thus we will 
take this value as the scale of coarse graining in the QCD 
dynamics (i.e., the mean- free path of partons) to obtain 
the hydrodynamics of QGP (h ~ 0.1/m). 

The second step is to approximate this coarse grained 
distribution acc?( r ) t) by replacing the integral in Eq. (Tl9|) 
by a summation over a finite and discrete set of points, 
{r a (t),a = 1,..,N SP h}, 



Nsph 

E 



A a (t)W(\r-v a (t)\). 



(21) 



If the choice of {A a (t), a = 1, .., Nsph} an d 
{r a (t), a = 1, .., Nsph } are appropriate, the above 
expression should converge to the coarse-grained 
distribution ace f° r large Nsph- Parameters 
{A a (t),a = l,..,NsPH} and {r a (t),a = l,..,NsPH} 
should be determined from the dynamics of the system. 
In practice, we first choose the reference density a* 
which is conserved, 



da* 
~~dt 



+ V-j = 0, 



(22) 



where j is the current associated with the density a* 
Then, we note that the following ansatzs, 

Nsph 

a SPH (r,t) = J2 »aW(\r-r a (t)\), 

a=l 

Nsph , /,\ 

i SPH (r,t) = ]T vJ^-W{\v-v a {t)\l 



a=l 



satisfies the equation, 



da 



SPH 



dt 



v • j S PH = o, 



where v a 's are constant. By using the normalization of 
W, Eq.®, we have 



,t)d 3 



SPH 



Nsph 
a=l 



Then we can interpret the quantity v a as the conserved 
quantity attached at the point r = r a (t). Therefore, 
the distribution a* SPH (r, t) is a sum of small piece-wise 
distribution, carrying the density, 

v a W{\v-r a [t)\). 

These pieces are referred to as " SPH-particles" . 

Using the above reference density and the extensive 
nature of A, we can write A a in Eq. ((2T|) as 



A a (t) 



a(r ai t) 
'a*(r a ,t) 



which represents the quantity A carried by the SPH par- 
ticle at the position r = r a (t). In fact, the total amount 
of A of the system at the instant t is given by 

Nsph 

A(t)= J2 

a=l 

In the previous works the entropy density is chosen as 
the reference density and the dynamics of the parameters 
{v a {t),a = 1, .., Nsph} are determined from the varia- 
tional principle from the action of ideal hydrodynamics 
(2^ | . Thus, the SPH particle coordinates and their time 
derivatives are considered as the variational parameters 
which optimize dynamically the action of the system. 

The entropy density is, however, conserved only for 
the motion of an ideal fluid. When we discuss the behav- 
ior of viscous fluids, where the entropy is not conserved, 
we cannot use the entropy density as the reference SPH 
density. Thus we introduce a new conserved quantity, the 
specific proper density a, which is defined by the flow of 
the fluid, 



3„ K) = 0, 



(23) 



and we will use it as the reference density for viscous 
fluids. Here, the four-velocity is defined in terms of 
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the local rest frame of the energy flow (Landau frame). 
The specific density is expressed in the SPH form as 

Nsph 

a*(r,t) = Y, v a W(\r-T a (t)\), 

a=l 

where a* = cru is the specific density in the laboratory 
frame and v a is the inverse of the specific volume of the 
SPH particle a, and is chosen as an arbitrary constant. 
Final results do not depend on this choice and we set 
v a = 1 for simplicity. As for the kernel W(r), we use the 
spline function [25j . 

Strictly speaking, this procedure is only possible pro- 
vided that the lines of flow in space defined by the veloc- 
ity field u M do not cross each other during the evolution 
in time. That is, if there appear turbulences or singular- 
ities in the flow lines, the above definition of Lagrange 
coordinates fails. However, if the size of h is consistently 
chosen as the size of coarse graining of the underlying 
microscopic theory, the flux line calculated using this h 
should not cross. 

Now we apply this method to the causal dissipative 
hydrodynamics in f + 1 dimension. We have to solve the 
evolution equation of the viscosity in the SPH scheme. 
To do so, we express the viscosity as 



It should be noted that there exist ambiguities within 
the resolution of the coarse-graining size h to express the 
equation of motion in the SPH form. However, in the 
ideal fluid, the SPH equation of motion can be derived 
by the variational method uniquely. Thus, we obtain the 
equation of motion by using the same SPH parametriza- 
tion to Ea.([2T|). 



d ( e a + p a + Il a i 
cr a - — I u, 

dTr 



Co 



Nsph 
0=1 



Pfj + Up p a + U a 



diWiK-rpW), 



(28) 



where the right hand side of Eq. (|28|) corresponds to the 
term di {p + n) written in terms of the SPH parametriza- 
tion. We remark that in the case of vanishing viscosity 
our result is reduced to the expression derived with vari- 
ational principle for ideal fluids. 

By separating the acceleration and force terms in 
Eq. (|28p . we obtain our final expression of the equation 
of motion for the SPH particles, 



N S i 



n = E ^W{\v-v a {t)\), 



(24) 



du a 



Time evolution of the term fl Q can be calculated as 



(Mr 



7a" 



dt 



'—(d^) a -—U a (25) 
tr tr 



where j a is the Lorentz factor of the a — th particle. At 
the same time, using the SPH expression for the entropy 
density s* in the observable frame, 



Nsph 



•'a(t)|), 



(26) 



where 



with 
A 
D 



M i:j = 7 (e + p + U)Sij + Aui 



j ■ 



f 

a* at 



-dj (p + n) + Bu 



3 ' 



n- — (e+p) 
as 

n 

tr 



n 

s H 

T 



TR 



(29) 
(30) 



(31) 
(32) 



and using Eq.([3]) we find 

d 



dt KcrJa ~ Ta* [l 



where s — s* /u is the proper entropy density. In the fol- 
lowing, we denote the quantity in the observable frame 
with the asterisk. In the above expressions, relaxation 
time Tr, viscosity coefficient £ and temperature T are 
functions of space and time, so that they should be eval- 
uated at the position of each particle a. 

Finally, we need to express the momentum conserva- 
tion equation by the SPH variables. We write the space 
component of Eq. |T]) in terms of the reference density, 



dr 



P 



n 



d l {p + Tl)=0. 



(27) 



These set of equations define the coarse-grained dynamics 
to represent the continuity equation for the energy and 
momentum tensor of a relativistic fluid, together with an 
irreversible mechanism which converts a part of collective 
kinetic energy (the motion of SPH particles) into internal 
heat of the fluid. 

For the sake of book-keeping, we summarize the prac- 
tical algorithm of calculating the SPH method to solve 
numerically the causal dissipative hydrodynamics. The 
dynamics is described by the following variables corre- 
sponding to the quantities attached to the each SPH par- 
ticle: 

^r a ,Uo,,y a , (J^J ,U a ; a = l,..,N S PH^ ■ 

At the initial time, their values are determined accord- 
ing to the initial condition. The entropy density profile 
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and the bulk viscosity are then obtained with the in- 
terpolations in Eq. (j2U)) and Eq. (pM|) , respectively. The 
energy density, pressure and temperature are calculated 
with the equation of state. The time evolution of these 
quantities are calculated by solving the equations derived 
previously 

du a _ M -ip 

= --(W a --n., 



dt 



1(1) 

dt W< 



tr 



TR 



Here the inverse matrix of M is calculated explicitly as 
1 



M- 1 = 



7(e + p + n 



■1 



7(e + p + n) ( 7 (e + p + II) + ( 7 2 - 1) A) 
and the four-divergence of the velocity is calculated as 
7a ( dcr* " 



1 _ gK q 

7a 



with 



da* 
dt 



J- W* - V- j. 



1 
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FIG. 3: The time evolutions of the proper entropy density for 
ideal fluid (dotted line) and the viscous fluid of a — 0.1 (solid 
line) . 



which in turn, 



IV. EXAMPLES 

A. Expansion to the vacuum and stationary 
boundary 

Let us consider first the Landau type initial condition. 
For the ideal case, this example has already been dis- 
cussed in and the SPH scheme works very well. When 
we introduce the viscosity, we found that the sharp dis- 
continuity in the boundary leads to undesirable instabil- 
ities. As we will discuss later, the origin of these insta- 
bilities is due to the presence of a space discontinuity. 
When we relax such a steep boundary by replacing the 
initial distribution consistent with the SPH scale used, 
the above instabilities disappears. In this subsection, we 
use basically the Landau type initial distribution with 
the temperature 200 MeV, distributed uniformly within 
the range x G [—1,1] fm- To relax the sharp boundaries, 
we add the surface thickness of lOh, where h = 0.01 fm. 
For simplicity, we take both vanishing initial velocity and 
viscosity, 

u(x,0) = 0, II(x,0) = 0. 

As we mentioned in the introduction, we use the equation 
of state that of massless ideal gas, 



P : 



(33) 



where C is a constant related to the Stephan-Boltzmann 
constant of massless 3 flavor quark-gluon gas. 

The temporal evolution of the density profile for the 
entropy at t = 5, 10 and 15 fm is shown in Fig[3l The 
solid line represents the results of the causal hydrody- 
namics with small viscosity, a = 0.1 and b = 6. For the 
sake of comparison, we show the time evolution of the 
same initial condition for the ideal fluid with the dashed 
lines. 

For the massless ideal fluid, we know that the prop- 
agation speed into the vacuum should be the speed of 
light. On the other hand, one can see that the propaga- 
tion speed into the vacuum for the viscous fluid is slower 
than that of the ideal fluid. This is because the viscosity 
acts as an attractive interaction during the expansion of 
the fluid. Thus it takes more time to achieve the speed 
of light in the causal dissipativc hydrodynamics. On the 
other hand, we can also observe that the propagation of 
rarefaction wave into the matter of the viscous fluid is 
faster than that of the ideal case. This is what we expect 
from Eq. (fT8|) . 

Interestingly enough, we observe that the behavior of 
the boundary of the viscous fluid seems to be a kind of 
a stationary wave. To understand this phenomenon, let 
us introduce the relative coordinate y = x — v s t with 
a velocity of the stationary wave v s and suppose that 
the energy-momentum tensor depends only on y near the 
boundary. Then, from the equation of continuity of the 
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360 



FIG. 4: The time evolution of the pressure P (dotted) and 
viscosity II (solid) of the fluid. At the boundary, we can see 
the relation P — —II. The in-set shows the details of the 
region where the two curves starts to coincide. 



energy-momentum tensor, we have 



d 

dy l 
d_ 

dy l 



— \-v s T 00 + T 01 } = 0, 



— \-v.T 01 +T 11 ] = 0. 



(34) 
(35) 




where, 



FIG. 5: The shock formation in the proper entropy profile 
(solid) of the viscous fluid of a = 0.1, starting from the dis- 
continuous initial condition (dotted). 



We have checked for various combinations of parame- 
ters and we found that the above feature of the appear- 
ance of a stationary wave near the boundary is universal 
for the free expansion of viscous matter into the vacuum. 
This fact is very important because the viscosity II is 
always the same order as the pressure P near the bound- 
ary to the vacuum, showing that the viscous effect can 
not be treated as a small correction to the equilibrium 
thcrmodynamical quantities. 



I* 4 " = (e + P + II) vfu" - gT (p + IL). 

and — (y, •yv) . From the boundary condition at the 
vacuum where vanishes, we get the following two 
conditions to define the stationary wave, 



v s (T 



li 



rp00\ 
Jill 



T 01 (l- 
v 2 e T 00 . 



(36) 
(37) 



From Eq. (|36p . we find v s = v or l/v. The first solution 
is consistent of our stationary wave assumption, showing 
that the fluid propagates to the vacuum with the veloc- 
ity v. The second solution leads to a trivial situation 
e = P = 0. When we substitute v s = v into Eq. ([3"T]) we 
get P = —II. Thus, if the boundary of the fluid prop- 
agates as a stationary wave, then the pressure and the 
viscosity should satisfy this relation. In fact, this relation 
corresponds to the case where the acceleration vanishes 
for the fluid motion. In FigHJ we plotted the profiles of 
P and II and we can see clearly the relation P = — H is 
satisfied near the boundary (see the in-set). In this fig- 
ure, we notice that there appear quick oscillation modes 
with the wavelength of the order of few h. We will discuss 
this point in the following sections. 



B. Density discontinuity and shock wave 

One interesting question is whether the dissipative hy- 
drodynamics can describe the formation of a shock wave. 
We know that for an ideal fluid, when a shock is formed, 
we need the so-called artificial viscosity for smoothing 
the shock region. One might argue that when the real 
viscosity is present, we do not need to introduce such an 
artificial viscosity. In Fig. [5l we show the time evolution 
of the entropy density of a viscous fluid of a = 0.1 whose 
initial distribution has a discontinuity. We know that for 
the case of an ideal fluid, such an initial distribution leads 
to a shock wave propagation into the medium. The thick- 
ness of the shock front should be zero for an ideal fluid. 
With the presence of a finite viscosity, the shock wave 
is smoothed and appears just a quick change in density. 
For the comparison, we show in Fig. [5] the time evolution 
of the density profile for an ideal case calculated by the 
SPH scheme with a finite h, where we see a quick oscilla- 
tion of large amplitude. It is interesting to observe that 
the wavelength of these oscillations is exactly the order 
of h of the SPH scheme. 

The physical reason for these oscillations of the SPH 
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FIG. 6: The behavior of SPH solution for the ideal fluid in 
the presence of shock wave corresponding to Fig. [5] 



FIG. 7: The shock formation in the proper entropy profile 
(solid) of the viscous fluid of a = 1, starting from the discon- 
tinuous initial condition (dotted). 



particles is clear. From the well-known Hugoniot- 
Ranking relation, we know that at the shock front, there 
should exist a production of entropy (l8j |. However, for an 
ideal fluid, the total entropy is conserved so that the SPH 
particles carry the extra kinetic energy corresponding to 
the entropy production at the shock. These extra oscilla- 
tions propagate with the smallest wavelength (order of h 
). The presence of viscosity can damp these oscillations 
and turns the kinetic energy of the SPH particles to in- 
ternal heat of the fluid, recovering the Hugoniot- Ranking 
relation. This is the case of the example in Fig. 

However, this is not always the case. For example, 
when a becomes large, the relaxation time tr increases, 
because of causality. Then the time scale for the damping 
becomes comparable to the evolution time scale generat- 
ing an oscillatory behavior. Such example is shown in 
Fig. [7] with o=l. In addition, when we look precisely 
the density profile shown in the example of Fig. we 
note the existence of quick oscillation modes with the 
wavelength of the order of h, although in this example, 
their amplitude does not increase in time. 

In the SPH calculation, by assumption, h is the scale of 
coarse graining and the resulting dynamics should always 
have larger wavelength than this. Thus, it is clear that 
the appearance of such rapidly oscillating modes of short 
wavelength of the order of h is a signature of an inconsis- 
tency of the theory with this finite coarse-graining scale. 
Such a situation also happens in a simple expansion of 
the fluid into the vacuum in the presence of viscosity. As 
mentioned in the beginning of this section, if we take a 
real sharp discontinuous Landau initial condition, sim- 
ilar rapidly oscillating modes appear at the expanding 
boundary. This is because a real discontinuity is not com- 
patible with a finite coarse-graining scale h. See also the 
small oscillations in Fig. 0J When we discuss the shock 



phenomena using an ideal fluid, then the shock front has 
null thickness and it should be treated as a singularity 
of the theory. This is because the usual hydrodynamics 
assume a null coarse-graining scale. 

Note that the usual viscosity is also obtained assuming 
a vanishing coarse-graining scale. To be consistent with 
the coarse-graining scheme of the theory, there should 
be an additional dissipation mechanism which is related 
to the coarse-graining scale. Such a mechanism should 
vanish in the limit of h/L <C 1, where L is the typical 
observable scale of the dynamics. 



V. VISCOSITY ASSOCIATED WITH COARSE 
GRAINING SCALE 

As discussed above, we need an additional viscosity, 
whose scale is determined by h. Such a viscosity also 
should obey the requirement of causality. Thus we pro- 
pose 

n Tot = n + m, 



where as before 



.(h) 



dtt h 
dr 



with 



-(h) 



-(h) 



(e + P) h, 



e + p 



The coefficients and are the numbers of the or- 
der of unity. Strictly speaking, these numbers should be 
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determined from a microscopic theory by incorporating 
the effect of the coarse-graining scale appropriately. Here 
we take them as phenomenological parameters and found 
that the values, 



2GC 



1/2, &W=2, 



(38) 



can eliminate undesirable oscillations in the SPH motion. 

When this additional viscosity is included, the viscosity 
term IT in Eqs. (f29|) and (f30|) should be replaced by Hrot 
and the coefficients A and B in Eqs. ([3T|) and (j32|) are 
modified as 



A 



e + p- 



n 



Tot 



d {e + p) 



ds 



s + 



n 



Tot 



T 



TR 



-(h) 



B = A 



7 2 da* 
a* dt 



n 

TR 



(39) 
(40) 



Consequently, the dispersion relation for the linearized 
sound wave, Eq. lfTTj) is modified as 



-ak 1 



C 



L0k 2 



+ 1- 



[e + p)l + vtrw (e + p) i + i T 



(h), 



(41) 



In Figs. [8] and El we show the real and imaginary 
parts, respectively, of the frequency u> as function of k, 
calculated from Ea.([3I]) for a = 0.1, b — 6 and T = 
200 MeV. This time, there are 4 solutions, but only 2 
of them are propagating modes and the other 2 are non- 
propagating modes. All them have positive imaginary 
part so that our equation is stable to a linear perturbation 
around a hydrostatic state. 

Under the normal condition, t^'/tr -C 1, "C Ci 
the above equation gives the same dispersion relation as 

the previous case, as far as since the second 

term of the right-hand side of Eq. (|4"Tj) can be neglected 
compared to the first term. On the other hand, for the 
very large k limit, the asymptotic form of the dispersion 
relation becomes 



±\ a + - 



1 



1 1/(67*) 

l 2 a + 1/6 



1/(6^ 



.(hh 



1/6W 



The most effective choice of the additional viscosity is 
obtained for the smallest possible value of rjj , preserving 
causality. This determines the value of b^ h > as 



1-a-l/b 

The choice Eq. (|38)) corresponds to the case 6 = 6 and 
a = 1/3. The difference of the two regimes determined 
by each viscosity can be seen in these figures. 

The additional viscosity introduced here can be consis- 
tent with the stability and causality of the theory. With 
this procedure, the rapid oscillation modes associated 
with the degrees of freedom beyond the applicability of 
the coarse-grained theory can be naturally eliminated. 




80 120 
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FIG. 8: The real part of the frequency ui as a function of k. 
There are two propagating modes A, B (dashed lines) and 
two non-propagating modes C and D which are degenerated 
(solid). 
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FIG. 9: The imaginary part of the frequency u> as a function 
of k. The two propagating modes A and B are degenrated 
(dashed line). The non- propagating modes C and D have 
different k dependence. 



A. "Double Shock" Phenomena 

We have already shown that, without the additional 
viscosity, the SPH time evolution develops rapidly oscil- 
lating modes with large amplitude whose wavelength are 
of the order of h. In Fig. [TQl we show the result of the 
same time evolution with the additional viscosity. We see 
that this new viscosity successfully smears out the quick 
oscillating modes. 
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FIG. 10: The same example as Fig. [7] calculated with the 
additional viscosity. Quick osculations at the shock region are 
removed. 



FIG. 12: The shock formation of the proper entropy density 
of the ideal fluid, starting from the homogeneous initial con- 
dition (dotted). 
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FIG. 11: The velocity profile corresponding to Fig. 1101 



A very interesting feature of this example, in Fig. 1101 is 
the appearance of non-trivial structure of discontinuities 
when compared to the usual shock discontinuity seen in 
Fig. [4j By changing the viscosity parameter a from the 
small value to the large ones, we found that the discon- 
tinuity associated with the usual shock corresponds to 
the discontinuity A shown in this figure. The new dis- 
continuity B starts to appear at some critical value of 
parameter a. This can be interpreted as the transporta- 
tion of the discontinuity in the initial distribution by a 
causal diffusion mechanism, as is known in the case of 
matter transportation by a telegraph equation (See Ap- 
pendix ). In fact, the velocity profile for this apparent 
"double shock" phenomena does not have any disconti- 
nuity at the location of B as is shown in Fig. [TT1 



B. Shock wave Formation by quickly moving fluid 
component 

Another typical example where shock discontinuities 
appear is that when a component of fluid is accelerated 
by an external force and achieves a velocity greater than 
that of the sound in the medium. Such phenomena are 
particularly interesting in heavy ion physics. When high 
energy partonic jets punch out the thermally equilibrated 
QGP, they may transfer a large amount of momentum 
and energy, dragging the piece of the fluid with high ve- 
locity. Such a scenario is often discussed with the ob- 
served angular correlation of produced hadrons, connect- 
ing to the formation of Mach cones[3,0|. Of course, the 
Mach cones do not exist in 1+1 dimensional systems, 
but it is interesting to study the dynamical formation of 
a shock discontinuity in our theory. If we calculate such 
situation without the additional viscosity, the formation 
of shock leads to quick unphysical oscillations as shown 
in Fig. [T2] Here, a small part of the fluid has an finite ini- 
tial velocity (Lorentz factor 7 ss 2) with constant entropy 
density in the local rest frame. However as is discussed 
in the previous section, the inclusion of the additional 
viscosity can eliminate these unphysical oscillations as 
shown in Fig. I13[ confirming the efficiency of our addi- 
tional viscosity. Differently from the case of discontinuity 
initial condition shown in Fig. [4l we found that the shock 
front created by a rapidly moving fluid element cannot 
be smoothed out even in the presence of normal viscos- 
ity. The effect of piling up the matter at the shock front 
generates a steep density variation whose wavelength be- 
comes eventually smaller than the coarse-graining scale 
h. Therefore, the introduction of the additional viscosity 
is essential to calculate the dynamical shock formation 
processes. 

In Figs. [FHand[T31 we show the entropy density pro- 
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FIG. 13: The same as Fig. 1121 calculated with the additional 
viscosity. 
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FIG. 16: The corresponding velocity profile of the viscous 
fluid of a — 1. 




FIG. 14: The same as Fig. [T3]with a = 0.1. 
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FIG. 15: The same as Fig. [TJwith a = 1. 



file of the same situation for different values of viscous 
coefficient a with additional viscosity. It is interesting to 
note that a rarefaction wave is formed and propagates 
backwards with a smaller velocity, as is shown in these 
figures, although the fluid velocity is positive everywhere, 
as is shown in Fig. [TBI Note that, another density dis- 
continuity is formed in the rarefaction wave part which 
behaves exactly in the same manner as the case of shock 
front created by the initial density discontinuity. 



VI. DISCUSSION AND CONCLUDING 
REMARKS 

In this paper, we studied the causal dissipative hydro- 
dynamics in 1+1 dimensional systems. To clarify the 
effects of viscosity and its relation to the coarse-graining 
aspect of the theory, we apply the SPH formulation to 
represent the macroscopic variables. Here, we consider 
the SPH scale parameter h as the scale of coarse-graining. 

We have shown that once the viscosity is introduced 
appropriately, then the expansion of the fluid into the 
vacuum should form a steady wave and there the univer- 
sal relation p + II = should be satisfied. We also stud- 
ied the various situations where the shock discontinuities 
emerge. We argue that for the consistency of the causal 
dissipative hydrodynamics, we have to introduce an addi- 
tional viscosity associated with the coarse-graining scale 
of the theory. This is because, the normal viscosity ob- 
tained from microscopic theories, such as kinetic equa- 
tions or the Green-Kubo-Nakano formula, is associated 
with the bulk properties of the matter and is calculated in 
the vanishing coarse-graining scale h — > 0. To be consis- 
tent with the coarse-graining scheme of the theory, there 
should be additional dissipation mechanisms that disap- 
pears in the limit of vanishing h. We proposed a scheme 
which contemplates such an additional viscosity, keeping 
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causality and stability of the theory. 

Most of the results in this study will be more relevant 
when the role of viscosity becomes effective. Such situ- 
ations will be realized in the coming LHC experiments 
where we may expect that a large inhomogeneity in the 
velocity profile can be created. The application of the 
present theory for the LHC energy is now in progress. 

One might think that the problems discussed here 
might be out of the range of applicability of hydrody- 
namics in the usual argument based on the smallness of 
the Knudsen number. However, as we already pointed 
out, it is well-known that hydrodynamics is still applied 
even for cases where the deviation from equilibrium is 
apparently large. In fact, the microscopic derivation 
of hydrodynamics is still an open question. Interesting 
works are still under progress to justify the applicabil- 
ity of hydrodynamics based on the asymptotic theory, 
the fluctuation theory and so on. In our opinion, while 
we do not establish a precise theoretical criteria for its 
use, hydrodynamics should be explored independent of 
such formal limitations, since it contains very important 
ingredients to describe phenomenologically the collective 
aspects of matter. As a matter of fact, if we had really to 
restrict ourselves to the region where the hydrodynamics 
approach is clearly available, we would never discuss the 
phase transition dynamics in this scheme [28| , so that all 
the ideal fluid model for heavy-ion collisions would be- 
come meaningless. In the same sense, we consider that 
it is significant to study the causal dissipative hydrody- 
namics comprehensively. 

In this paper, we discussed the problem of causality 
and stability of our theory, and concluded that the theory 
is stable for the small linear perturbations. To be precise, 
we did not prove that the theory is stable for nonlinear 
perturbations, although the numerical solutions exam- 
ined here did not show any instabilities. However, if the 
theory is shown to be unstable in some regime, then we 
have to construct another stable theory. Analysis on this 
line is under investigation and will be reported in another 
paper. The effect of the equation of state containing 



phase transitions also will be discussed in future. 
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APPENDIX A: TELEGRAPH EQUATION 

We discuss the behavior of the causal diffusion equa- 
tion, which is described by 



a 2 



d 



n d2 



0. 



As is well-known, the maximum propagation speed of 
signals described by the equation is v — y/D/r. As is 
discussed in [29( , the solution of the causal diffusion equa- 
tion is given by 



n(x, t) — —e 2t [no(x + vt) + no(x — vt)] 




x^-n (z ), 
at 



where Uq(x) is the initial distribution and Iq(x) is the 
modified Bessel function. The first term of the solution 
indicates that the initial distribution iiq(x) is separated 
into two fragments, which travel to the opposite direc- 
tions with the velocity v. Thus, in the shorter time scale 
than the relaxation time r, the memory of the initial 
distribution profile survives near the boundary. 
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